function [inten_full] = intensity(p,dp,q,dq,V,alpha,sigma,gamma)
    
    xi = 2*alpha/sigma;
    I = size(V,1);
    J = size(V,2);

    % Vectorize
    p = p(:);
    q = q(:);
    V = V(:);
        
    % Estimate Vpq (includes reflecting boundary conditions)
    injndiag = ones(I*J,1)/(4*dp*dq);
    ijndiag = zeros(I*J,1);
    ipjndiag = -ones(I*J,1)/(4*dp*dq);
    injdiag = zeros(I*J,1);
    ijdiag = zeros(I*J,1);
    ipjdiag = zeros(I*J,1);
    injpdiag = -ones(I*J,1)/(4*dp*dq);
    ijpdiag = zeros(I*J,1);
    ipjpdiag = ones(I*J,1)/(4*dp*dq);
    
    A = sparse_matrix_function(injndiag,ijndiag,ipjndiag,injdiag,ijdiag,...
        ipjdiag,injpdiag,ijpdiag,ipjpdiag,I,J);
    
    Vpq = A*V;
    
    % Estimate Vpp (includes reflecting boundary conditions)
    injndiag = zeros(I*J,1);
    ijndiag = zeros(I*J,1);
    ipjndiag = zeros(I*J,1);
    injdiag = ones(I*J,1)/(dp^2);
    ijdiag = -2*ones(I*J,1)/(dp^2);
    ipjdiag = ones(I*J,1)/(dp^2);
    injpdiag = zeros(I*J,1);
    ijpdiag = zeros(I*J,1);
    ipjpdiag = zeros(I*J,1);
    
    A = sparse_matrix_function(injndiag,ijndiag,ipjndiag,injdiag,ijdiag,...
        ipjdiag,injpdiag,ijpdiag,ipjpdiag,I,J);    
    
    Vpp = A*V;
    
    % Grid search across state space
    n = 0:0.0001:1; % intensity grid
    cdim = size(n,2);
    rdim = size(Vpp,1);
    
    n_mat = repmat(n,rdim,1); % search for intensity in columnular dimension
    Vpp_mat = repmat(Vpp,1,cdim);
    Vpq_mat = repmat(Vpq,1,cdim);
    p_mat = repmat(p,1,cdim);
    q_mat = repmat(q,1,cdim);
    
    f = -gamma*n_mat.^(gamma-1) + 0.5*Vpp_mat.*(p_mat.^2).*(1-p_mat.^2)*xi^2 + ...
        0.25*(Vpq_mat.*p_mat.*(1-p_mat).*q_mat.*(1-q_mat).*xi^2)./sqrt(1+n_mat);
    [~,ind] = min(f.^2,[],2);
    
    ind = sub2ind(size(f.^2),(1:size(f.^2,1))',ind); % need to convert to linear index
    inten = n_mat(ind);
            
    inten(p==0) = 0; % impose zero info acq at boundaries
    %inten(p==1) = 0;
    inten(q==1) = 0;
    
    ind = (p>0.98);
    inten(ind) = 0;

        
    % reshape to IxJ matrix
    inten_full = reshape(inten,I,J);
 
    